This document describes fits from both state-space Ricker spawner-recruit models used fit in the Research Document (link once published online):
citation
Two state-space spawner-recruit models were fit in this paper that serve different purposes. First, a Ricker model with autoregressive (AR1) recruitment residuals was fit to estimate biological benchmarks (\(S_{gen}\), \(S_{MSY}\), and \(U_{MSY}\)), then a model with time-varying productivity (Ricker \(\alpha\) parameter) was fit in order to condition the biological submodel on recent population dynamics in the closed-loop forward simulation. These models were fit to each of the 9 CUs one by one, yielding 18 models to diagnose total. Data and code to reproduce the analysis and report is available in this GitHub repository, where you can see the R code to run the models, and the models themselves.
We fit the spawner-recruitment model in a Bayesian estimation framework with Stan [@carpenter_stan_2017; @standevelopmentteamRstanInterfaceStan2023], which implements the No-U-Turn Hamiltonian Markov chain Monte Carlo algorithm [@hoffman2014] for Bayesian statistical inference to generate a joint posterior probability distribution of all unknowns in the model. We sampled from 4 chains with 2,000 iterations each and discarded the first half as warm-up. We assessed chain convergence visually via trace plots and by ensuring that \(\hat{R}\) (potential scale reduction factor; @vehtari2021rank) was less than 1.01 and that the effective sample size was greater than 400. Posterior predictive checks were used to make sure the model returned known values, by simulating new datasets and checking how similar they were to our observed data.
These should be clearly mixed, with no single distribution deviating from others (left column), and no chains exploring a strange space for a few iterations (right column).
We hope minimum effective sample sizes are greater than 2000 and that \(\hat{R}\) are less than 1.05.
For the AR1 model:
| CU | ESS | Rhat |
|---|---|---|
| NorthernYukonR.andtribs. | 1022 | 1.002 |
| Whiteandtribs. | 525 | 1.010 |
| Pelly | 603 | 1.004 |
| Stewart | 1212 | 1.002 |
| Nordenskiold | 1674 | 1.003 |
| YukonR.Teslinheadwaters | 1110 | 1.003 |
| MiddleYukonR.andtribs. | 969 | 1.004 |
| UpperYukonR. | 1188 | 1.001 |
| Big.Salmon | 698 | 1.006 |
and the time varying model:
| CU | ESS | Rhat |
|---|---|---|
| NorthernYukonR.andtribs. | 70 | 1.050 |
| Whiteandtribs. | 572 | 1.008 |
| Pelly | 507 | 1.005 |
| Stewart | 48 | 1.091 |
| Nordenskiold | 1112 | 1.007 |
| YukonR.Teslinheadwaters | 813 | 1.003 |
| MiddleYukonR.andtribs. | 637 | 1.012 |
| UpperYukonR. | 62 | 1.056 |
| Big.Salmon | 465 | 1.016 |
We see some problems with the time varying model, and have explored various parametrizations and options to include semi-informative beta priors that did not improve diagnostics. We imagine these issues in the time-varying model are due to the extreme, pronounced decline in productivity, and the correlation between productivity and the stationary capacity prior (i.e., Ricker \(\beta\)).
**add a head(max(Rhat)) type thing to show what model*parms had problems?**